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Abstract 

The 3PI method is a technique to resum an infinite class of diagrams, which may be useful 
in studying nonperturbative thermodynamics and dynamics in quantum field theory. But it has 
<^ never been successfully applied to gauge theories, where there are serious questions about gauge 

invariance breaking. Wc show how to perform the 3PI resummation of QCD in 3 Euclidean 
spacetime dimensions, a warmup problem to the 4 or 3+1 dimensional case. We present the 
complete details of the technical problems and how they are overcome. We postpone a comparison 
^ of gauge invariant correlation functions with their lattice-determined counterparts to a future 

publication. 
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I. INTRODUCTION 



The early Universe existed in a state of deconfined quark-gluon plasma, a state which has 
also recently been produced in the laboratory via heavy ion collisions. The thermodynamics 
of such plasmas are now relatively well understood. In the early Universe the gross features 
are well described by perturbation theory, and while the strength of the electroweak phase 
transition (or crossover) cannot in general be determined perturbatively, there are powerful 
lattice methods which can be brought to bear p]. Lattice methods can also describe the 
thermodynamics of QCD near its intrinsic scale Aqcd, provided that the chemical potential 
is not too large |2H1]. However, many of the most interesting questions in each context have 
to do with dynamics or unequal time correlations. For instance, real-time correlations and 
nonequilibrium currents are important in understanding whether electroweak baryogenesis 
can occur at a first order electroweak phase transition [51 E]. In the heavy ion context 
there are many dynamical quantities we would like to know, such as the viscosity [7HTU]. 
heavy quark diffusion rate [HI [12], photon production rate [T314T5] . and so forth. Dynamical 
properties of QCD are also important at much higher temperatures such as the electroweak 
temperature, where they could play a role in baryogenesis and in various phase transitions. 

The problem is that we have no first-principles, intrinsically nonperturbative technique 
for theoretically predicting such real-time properties, even in equilibrium or linear response. 
We have models and phenomenological fitting (for instance, of the viscosity using elliptic 
flow in heavy ion collisions PUj), but the only tool we have which is close to first principles 
is perturbation theory. 

Perturbation theory is notoriously poorly convergent when applied to hot gauge theories. 
For instance, the expansion for the pressure of QCD as a series in as is known to order 

In as, but the series only appears to be useful at temperatures many times the scale Aqcd 
[16] . For real-time quantities, certain (Hard Thermal Loop or HTL) resummations are 
necessary even to find leading order results for transport coefficients [18]. Even so, the per- 
turbative expansion for real-time quantities appears to be even worse behaved than it is for 
the pressure and other thermodynamic quantities, at least if we restrict attention to trans- 
port coefficients, which involve either zero frequency and momentum or lightlike momentum 
limits of external 4-momenta. Only two such quantities are known beyond leading order; 
the diffusion coefficient for a heavy quark [19] and the transverse momentum diffusion for 
a fast charge [20] ■ In each case the first corrections enter at order g, not order g"^ as would 
be normal in a vacuum perturbative expansion, and the corrections represent of order 100% 
shifts in the transport coefficients for as ~ 0.05, a value obtained only at temperatures well 
above 1 TeV! Though it is a little dangerous to extrapolate from two examples, it appears 
that even at the electroweak scale, the QCD sector of the Standard Model is probably not 
well described by perturbation theory, as far as dynamics are concerned. 

Let us analyze the problem in a little more detail, for the case of heavy quark diffusion 
[T^ . The leading-order diagram and the four next-to- leading order diagrams are shown 
in Fig. [T] Here the double line is the Wilson line heavy quark trajectory, and hatched 
blobs represent Hard Thermal Loop (HTL) resummation. In particular, the leading-order 
diagram is already one-loop self-energy resummed. The contribution from diagram A is 5 to 
10 times as large as that of the other diagrams (in Coulomb gauge, the choice used in [T9]). 
Therefore the physics which is problematic for the perturbative expansion is presumably 
the physics represented by this diagram. This diagram resembles the leading-order diagram, 
except that in the leading diagram the momentum in the self-energy loop is assumed to 



2 




A B C D 

FIG. 1: NLO heavy quark diffusion. 




FIG. 2: Inclusion of additional HTLs. 



be large compared to the propagator momentum, while in diagram A it is allowed to be 
of the same order. When the momentum is soft, it is also necessary to include HTLs on 
the propagators and vertices of the self-energy in diagram A. The large size of diagram A 
indicates that the soft loop momentum region in the self-energy is almost as important as the 
region where the loop momentum is hard. This is a problem because in this soft region, the 
HTL corrections must also be included on the propagators and vertices in diagram A. But 
why shouldn't these loops also receive large corrections from their soft-momentum regions? 
Then we also need to include the region where we replace the HTLs in diagram A with soft 
momenta, as suggested in Fig. |2} But when these momenta are soft, there are new vertices 
and propagators requiring HTL corrections. These may also have large soft corrections, 
bringing in more diagrams - and we are "off to the races." 

The problem is that the low-momentum {gT) region is not really a small part of phase 
space, and it is not really weakly coupled. But the good news is that the dominance of 
diagram A above suggests that it is really only the soft corrections to the HTLs present 
in the previous order which are important. This leads to an infinite number of diagrams 
contributing, but only a restricted (infinite) set of diagrams being "most important" - those 
suggested at in Fig. [2j If these diagrams could be resummed somehow, then we would likely 
capture all the most important corrections, and the range of validity of perturbation theory 
might be significantly expanded. 

There is a procedure for performing an iterative resummation of all one-loop self-energy 
and 3-point vertex corrections, as suggested in Fig. [2j It is the 3-loop, 3-particle irreducible 
(3PI) resummation scheme |2I1[22]- The above discussion suggests that a 3PI resummation 
may capture the most important higher-order physics and greatly improve the convergence 
of the perturbative expansion - even capturing some nonperturbative information. But this 
is by no means guaranteed, since the 3PI approach performs an incomplete resummation of 
higher-order effects, which furthermore is not gauge invariant. We would like, eventually, to 
perform a 3-loop 3PI resummation of real-time, 3-1-1 dimensional QCD (or the full Standard 
Model) at small or intermediate coupling. However, as yet the 3PI technique has never been 
applied to nonabehan gauge theories. And we would also like to test whether the resum- 
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mation technique is effective and reliable in a somewhat more controlled setting. Therefore 
we feel it is necessary to consider a slightly simpler problem as a warm-up exercise, and 
as a testing ground for whether the 3PI technique is effective in nonabelian gauge theory. 
Therefore, in this work we will investigate the 3-loop, 3PI resummation of 3- dimensional 
(3D), SU(3) Yang-Mills theory, an endeavor we will call the "all threes" problem (3D, 3-loop 
3PI for SU(3)). 

Working in 3 Euclidean dimensions simplifies our problem in two ways. First, the UV 
behavior of 3D Yang-Mills theory is much milder than in 4D, since the theory is super- 
renormalizable. The second simplification is that the vacuum Euclidean theory has fewer 
Lorentz invariants than the finite-temperature, Minkowski theory (or the Euclidean theory 
with periodic time direction). For instance, in (vacuum) Euclidean space the propagator 
G^'^ is built from two tensorial structures and is a function of one invariant; G^'^{p) = 
T^'^Gt{p'^) + Ij^'^ G lIp"^) . In real time at finite temperature, or in imaginary time with 
periodic boundary conditions, it has more tensorial structures and is a function of two 
variables, p^ and the energy p^. 

Note however that 3D Yang-Mills theory is far from trivial. The flip side of super- 
renormalizability is that, since the gauge coupling is dimensionful, it establishes a scale 
(momentum scale p ~ (7^ or length scale A ~ where we expect strongly coupled, 

nonperturbative behavior. Therefore 3D Yang-Mills theory displays both weak or strong 
coupling, depending on the energy scale. 3D Yang-Mills theory is also physically interesting. 
At the thermodynamic level, the infrared behavior of thermal 4D Yang-Mills theory (with 
any fermionic matter content) is 3D Yang-Mills with an adjoint scalar a slight 

extension of the problem we consider. The nonperturbative scale = q^y) corresponds to 
the scale gf^T of the full theory (at leading perturbative order). It is believed that the 
poor convergence of perturbation theory at intermediate couplings in QCD is due to the 
nonperturbative physics of the 3D theory. Therefore in a sense studying the 3D gauge 
theory by nonperturbative means is treating most of the physics which makes thermal QCD 
poorly behaved (at intermediate couplings). 

The goal of this paper is to produce a complete solution to the 3 loop truncation of the 
3PI effective action for 3D QCD. As outlined above, this should be viewed as a warm-up 
problem to what we would really like to do, which is to apply the resummation in 3+1 
dimensions in a thermal (or even nonequilibrium) context. We feel that this first step is well 
motivated. As already emphasized, the 3D theory is a subset of the 3+ ID theory (the theory 
we really want to solve). And the 3D theory can be studied nonperturbatively on a lattice, 
which means we will be able to test the 3PI resummation procedure in a nonperturbative 
context by seeing whether its predictions are successful. Such a test is necessary because the 
3PI resummation captures only an incomplete and gauge-noninvariant subset of diagrams, 
so there is no hard guarantee that it will successfully reproduce the nonperturbative IR 
physics. This paper will concentrate on the resummation; we will return to the comparison 
with nonperturbative lattice studies in a follow-up publication. 

II. THE 3PI EFFECTIVE ACTION 

The 3PI effective action is obtained by a Legendre transform of the generating functional 
of connected diagrams 

W[J,K,L] = - log / V<p e-(^M+-^^+l^<^'+i^<^-') (2.1) 
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with one, two and three-particle sources, J, K and L (our notation is schematic, with fields 
generically denoted by and matrix indices and space integrations suppressed). Conjugate 
to J, K and L are the variables 0, G and G3, which are labelled as such since it can be proven 
that they are equal to the connected one, two and three-point functions. The functional 
derivatives of W are given by 

^ = ^(G3 + 3G0 + 0^), (2.4) 
and the 3PI effective action follows, 

5W 5W 5W 

r[0, G, G3] = j-^+ k'j^ +^ir~ ^t"^' ^' 



From Eqs. (2.3) and (2.4), we see that W generates connected as well as disconnected 
diagrams^, so it does not, strictly speaking, generate cumulants. However, W is still equal 
to the logarithm of a moment generating functional, hence once can expect that it as well as 
its Legendre transform F have well defined extrema. In terms of F the equations of motion 
for (j), G and G3 read 

£ = ^ = ^ = 0. (2.5) 

6(j) 5G 5G-i 

We will now specialize this procedure to the case of QCD. Since we consider pure-glue 
QCD there is no Higgs mechanism and we expect the one-point functions to vanish at the 
extremum^; hence we can set the VEVs of A^, c and c to zero and work only with the 
two and three-point functions. To write down the exact form of F relevant to the QCD3 
field content, it is useful to define the propagators and vertex functions^ listed in Table |l| 
The one-particle-irreducible 3-vertex V is related to the three-point function via G3 = G^V; 
extremizing with respect to G and G3 is equivalent to extremizing with respect to G and V, 
which we find to be more convenient variables. The 3PI effective action is given in terms of 
G, A, V and V by 



^ For instance, ST/SK generates (j)^ = {(j))"^ ■ 

^ When Grassman fields such as ghosts are present, the extremum is generally a saddle-point rather than a 
maximum. 

Throughout this work, we will generally suppress color indices and will frequently also suppress Lorentz 
indices, when the indices can be inferred from context. For instance, G — Gf^i, — GJ^^, VJ^i^p — V^^^p etc. 
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Resummed 



Bare 



Gluon Propagator 




Ghost Propagator 




Gluon 3- Vertex 


•frabc 


Ghost-Gluon Vertex 


\rabc 


Gluon 4- Vertex 





or ^^fr^ 



G 



{0)ab 



A(o)'''' = 



^^(0)abc 



2^M)abcd 




^fi^ o: 

or ^1 




or 




TABLE I: Feynman Rules for QCD3 



r = SiA^", c, c] + ^Tr log G-^Tt[G'^%^G-Tt log A + Tt[A'^%^ A 



+ 



+ 





e 




(2.6) 



(where we have explicitly written symmetry factors associated with diagrams and signs 
associated with ghost loops for clarity, as we will throughout), from which we have the 
following four equations of motion 



ST _ 6T _ 6T _ 6T 
SG ~ 6A~ W ~ 6Y 



0. 



(2.7) 



To illustrate the physics of Eq. (2.7), consider 5T/5G{p). Performing the variation using 
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the expression in Eq. (2.6), we find 



n(^)(p) 



(2.8) 
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which we recognize as the resummed one-loop self-energy (n = n(i) + n(2), we have not 
shown the similar graphical representation of 11*^^^). Similarly, variation with respect to V 



gives 






(2.9) 



which we recognize as a Schwinger-Dyson equation for vertex resummation. 

Note that the propagators appearing in all diagrams in Eq. ( 2.8[ ) and Eq. (2.9), as well 
as all vertices (except certain vertices in the one-loop self-energies), are the full objects. 
Therefore these equations must be solved self-consistently. The self-consistent solution of 
these equations represents our main challenge. We face two chief difficulties: 

• Decomposing the propagator into its transverse and longitudinal parts. 



Gti'{p)T^''' + (p)L 



^1/ 



p^p" 

p2 



g- 



(2.10) 



the propagators are determined in terms of three arbitrary functions of one continuous 
variable, Gt{p), Gl{p) and A(p) with p = \[p^. Similarly, the vertex V^i/i2A'3(Pi)P2,P3) 
(with = —p\y, — p2^i) can be expressed in terms of six independent tensorial struc- 
tures (see below), each multiplying an undetermined function of the three invariants 
pI, P2, and pi ■ p2 (or equivalently pf, P2, and pi). 

The challenge is that we are not merely solving for a few numbers, but self-consistently 
solving for unknown functions of one to three continuous variables. 

The one-loop gluon self-energy diagrams are linearly divergent, and the two-loop gluon 
self-energy diagrams are individually logarithmically divergent. These divergences 
must be regulated in a manner which respects gauge invariance,^ such as dimensional 
regularization. However since the propagators and vertices appearing in the diagrams 
are general functions of momentum which are presumably only known numerically and 
only in D = 3 dimensions, we will have to perform these integrations numerically. 



^ One might argue that, since the 3PI technique truncated to 3 loops is not gauge invariant, the use of a 
gauge invariant regulator is unnecessary. But we believe that it is necessary; first, our approach at least 
retains gauge invariance to low loop order, which would be lost without a gauge invariant regularization. 
And second, a gauge non-invariant regularization at 1-loop order would allow divergent masses, which 
fundamentally damage the physics. 
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The issue of divergences in self-energies is a technical issue which can be handled rather 
easily in 3 dimensions. The key fact is that at large momenta G,A,V and V approach their 
free values up to power suppressed corrections. We therefore know the exact form of the 
UV divergences. If we can find an expression with the same UV divergent behavior which 
is simple enough to integrate using dimensional regularization, we can add and subtract it. 
The subtraction renders the numerical integration of the full self-energy expressions finite, 
while the added version is integrated using dimensional regularization. We will explain this 
procedure in more detail in the next section. 

In order to fit arbitrary functions of one or a few real variables, we will write down a 
sufficiently flexible Ansatz for each function, with some set of variational parameters. That 
is, we take Gt{p) = Gxici^p) where q are coefficients - in practice, we take Gt{p) to be a 
rational function of p, and the q are coefficients of this rational function. Extremization of F 
with respect to G{p) is then replaced by its extremization with respect to the coefficients Cj 
of each propagator and vertex function. Under this procedure, some set of integral moments 
of the Schwinger-Dyson equations will be satisfied, rather than the equations being satisfied 
at every momentum value. We can determine the quality and limitations of this approach 
by seeing how the determined correlation functions change as the sizes of the variational 
Ansatze^ are changed; and we can directly test how well the Schwinger-Dyson equations 
are obeyed by computing directly the self-energies and vertex corrections at various mo- 
menta and comparing to the Ansatz, or measuring an integrated mean squared failure of 
the Schwinger-Dyson equations. 

III. DIVERGENCES AND REGULARIZATION 

As discussed above, variation of F with respect to a propagator gives rise to a Schwinger- 
Dyson equation involving self-energies written in terms of G and V. These will be, in 
general, complicated functions. Yet the self-energies may be UV divergent and so they must 
be regulated. Even after taming these self-energy divergences, the variation with respect 
to a propagator or vertex Ansatz coefficient SF/Sci may lead to a divergent integral over 
the propagator momentum p. We must also ensure that such divergences do not occur. We 
will handle these two problems in turn. Throughout we denote the momentum entering a 
self-energy as p, and use k and q for internal loop momenta. 

A. Divergences in self-energies 

In 3 dimensions the only divergent subdiagrams are gluon self-energies. To handle these 

divergences we must work in a regularization scheme which renders the self-energy diagrams 
finite and preserves gauge invariance. Therefore we will perform all integrals in dimensional 
regularization (DR), so J d^q J d^q with D = 3 -t- 2e. Unfortunately the self-energies 
contain the functions G and V, which are complicated and are only known in 3 dimensions. 
However, for any integral which is finite and well behaved in 3 dimensions, the e ^ limit 
of the DR value is the same as the value directly computed in 3 dimensions. Therefore we 
will start with identifying the UV divergent behavior of the full integrals containing G and 



^ The plural of Ansatz is Ansatze. 
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V so that we can subtract and add simple integrals with the same divergences. We can 
then perform the (finite) subtracted versions numerically in 3 dimensions, and finish off by 
adding back the simple integrals using DR. 

Dressed vertices and propagators are well behaved in the IR, hence the only divergences 
that we expect to see arise from the region of momentum space where q is large. Therefore we 
need to determine the asymptotic behaviors of G and V. Our theory is super-renormalizable, 
meaning that the coupling g"^ carries dimension, [g"^] = 1 = [g]. For large q, g"^ is small com- 
pared to the relevant scale q, so the large q region is weakly coupled and has a perturbative 
expansion. Further, powers of g^ in the expansion must be balanced against powers of q on 
dimensional grounds. Therefore the leading and first subleading behavior of the propagator 
in 3 dimensions is 

G,,{q) = ^(T,,(g) + eV(g)) + ^-^^T^^iq) + 0{q-') . (3.1) 
Similarly, the vertex goes as 

Vr^q + gy + 0{q-'). (3.2) 



The specific form of the g^ correction to V is known, see Appendix A3; but as we see in 
a moment we do not need it here. The one-loop correction to the gluon self-energy is also 
known. It is purely transverse and equals 



Ki'^iQ) = Q^-^ie + 2^ + ll)T,.(g), (3.3) 



hence 11^*^^) introduced in Eq. (3.1) is 



n^^^) = ^(e' + 2^ + 11). (3.4) 

64 



Now, consider the diagram 




(Pq 
(27r)3 



V,^sV,p.G''''ip + q)G'^iq), 



(3.5) 



where traces over internal color indices are implied (and hence the overall diagram is pro- 
portional to the color identity). Expanding the integrand in powers of q, for D = 3 the large 
q region of the integral behaves as 




2^2 



+ 2g\qr^\ + 2g\q)y^+ 0{g'q-') 



q'^ cp- 



+ ... 



NLO propagator 



V 

NLO vertex 



NLO^+NNLO 



(3.6) 

The first term arises from the leading order (bare) terms in the vertices and propagators, and 
the next two terms originate from the one-loop corrections (as marked). These first three 
integrals diverge, so we will have to add and subtract something to cancel their divergent 
behavior. 
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Actually, the NLO vertex corrections above will cancel when we sum over the one-loop 
self-energy corrections. To see this, consider the two diagrams, with one and with two full 
vertices: 






(3.7) 



The diagram with one full vertex enters with —2 times the weight of the diagram with two 
full vertices. Therefore the NLO vertex contributions from these two diagrams cancel, and 
the UV behavior is the same at NLO as the behavior of a loop with no vertex corrections. 
Provided that we perform the two diagrams by adding their integrands inside the integration, 
this cancellation takes place at the level of the integrand and does not lead to a log divergence 
in the integral in 3 dimensions. (This cancellation does not mean that the NLO vertex 
correction disappears; instead this correction will be accounted for explicitly when we include 
two-loop self-energy corrections.) 



Next consider the bare part of Eq. (3.6), which is linearly divergent in 3 dimensions. We 



will add and subtract a diagram made out of the bare vertex and propagator functions. 




K^SKSG(°)"^(p + g)G(°)^''(g) 



9 / y^as^vp 



(3.8) 



with G^^'^^s are y*^°)'s denoting bare propagators and vertices. The difference 





(3.9) 



is only logarithmically divergent in 3 dimensions. Moreover, Eq. (3.8) is finite when com- 
puted in DR and its D — )■ 3 limit is 




V- 



DR 
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(^(e' + 2e + ll)T^.(p)-(7^. 



p2 



(3.10) 



This, together with a contribution from the bare ghost loop which cancels the non-transverse 
piece above, gives rise to H^^^^ stated earlier. 

Lastly, we must subtract something with the same NLO "propagator" behavior remaining 



in Eq. (3.6). Naively, we could do this by defining 



GiY\p) = GS(p)n^(^'^)"'^(p)GS(p) = ^-5^T,.(p) 



(3.11) 



with Il^j^'''\p) and H^^^''') defined by Eqs. (A8) and (A9), and then "adding and subtracting" 
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the following diagram: 



A 



B 



p2 



+ finite . 



(3.12) 



The problem is that, as the above equation shows, the diagram is not only UV divergent (as 
expected, with coefficient A which we give below), but also IR divergent (with coefficient B, 
whose exact value will not be relevant). By including this diagram, we would be introducing 
an IR divergence where there is not one already. Rather, we will define an appropriately 
regulated self-energy corrected propagator. 



/i 



2e 



(p2 + m^) 2^ 



Pt^pu 



p2 _j_ ^2 



(3.13) 



so that the integration in 



e)IR Reg 5k 



(3.14) 



can still be performed analytically in DR. Upon integration, this diagram has the following 
form: 




A 



g^y + finite. 



(3.15) 



DR 



The coefficient that multiplies the UV 1/e, 



A 



7687r2 



(e + 4)(e2 + 2^ + 11), 



is identical to that of Eq. (3.12) due to the simple fact that 



hm Gi^/)^^^''^(g)= hm G^^^\q). 



(3.16) 



(3.17) 



Naturally, the finite parts of Eq. (3.12) and Eq. (3.14) will differ. 



The inclusion of such a gauge-noninvariant IR regulator is allowed because we are adding 
and subtracting this diagram, so any final results of our computation would be independent 



of the details of how this IR regulation is achieved. Returning to Eq. (3.6), we now have 



+ 





A 



~ —giMu + finite , 



(3.18) 



11 



The first line is the only part which contains full propagators; but it is finite at D = 3. 
Therefore its value in DR in the D — )• 3 limit simply equals its finite value in 3 dimensions, 
which we can find by numerical integration. The second line is divergent in 3D but can be 
carried out relatively easily in DR. It gives rise to the Ag^^je contribution, and to some 
of the finite terms presented in Eq. (A15) (Eq. (A15) is the sum of Eq. (3.14) as well as a 
similarly IR regulated Snowcone diagram). 

The procedure for handling the 2-loop graphs is similar. While the graphs are more 
complicated, the procedure is simpler, since in every case the only UV divergences arise 
when all components of the graph take their bare values. Further, no 2-loop graph we need, 
when built out of bare quantities, is IR divergent. Therefore we may simply subtract from 
each 2-loop graph, built using G and V , the same graph built using Go and Vq. The bare 
2-loop graphs can be performed in DR and will also give rise to a g^u/e divergence plus a 
finite part; the sum of these diagrams is given by Eq. (All). In the end, we find that all 
of the 1/e's cancel between two-loop diagrams: the 0{g'^) IR regulated gluon self-energy is 
UV finite in DR (the complete expression for which is given by Eqs. (All) and (A15)). 
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With this procedure in place, the Schwinger-Dyson equation, Eq. (2.8), becomes 



(3.19) 



+ 





1 

6 



FS*5> + 2 ■''SS'ii 



('{> 



n^(l'°)/ \ 1- /-l-rB(2,e)UV -p-rB(2,e)IR Reg 

where the n^^"'^) 's represent the sum all bare diagrams computed analytically in dimensional 
regularization (their values can be read from the Appendix, noting that Il^!)^''^^^^^°^{p) is 



exactly as stated in Eq. (A15), with no large p limit taken). In writing this equation we have 



suppressed the Lorentz indices, as written one should take 1X^^11'^'^ to get the transverse 
part of the self-energy needed to resum Gt and L^^^II^'^ to get the self-energy needed to 
resum Gl- 

All ghost self-energies are finite after angular averaging, and the vertex loops are power- 
counting finite, so no similar subtractions are needed in these cases. Nevertheless, we will 
still encounter divergences when it comes time to integrate these Schwinger-Dyson equations 
over propagator or vertex momenta. 
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B. Divergences upon p-integration 



Our plan is to find an approximate extremum of F by writing variational Ansatze for the 
propagators and vertices and to vary with respect to the Ansatz parameters. For instance, 
one could assume that the transverse propagator Gt{p) is the sum of a set of test functions 
with unknown coefficients. Grip) = Tli^i'Piip)- More generally, we choose Gt{p) to have 
some functional form with a set of variational parameters we will give our specific choice 



in Section IV Then variation of F with respect to Cj would yield 



and similarly for Gl and A. In the last section we ensured that the integrals involved in 
all self-energies 11 are finite. But this does not guarantee that the p integral above will be 
finite. For instance, if we chose the Ansatz 

Grip, example] = + j^^f^ + (^^^ (3.21) 

then SGt(p)/Sci = 2+^2 ■ And if ci ^ 1, then for large p, G^^{p) - G^^^~^{p) + IIt{p) ~ p^. 



In this case Eq. (3.20) would be cubically divergent. Physically this means that if we allow 
G{p) to vary from its correct value in a way which does not die away fast in the UV, then 
F will be divergently far from its extremum. 

Continuing with the same example, if we fix ci = 1, forcing the propagator to have the 
correct free-theory limit in the UV, then the most severe UV divergence arises from C2. For 
large p we have 6Gt{p)/Sc2 ~ and G^^{p) — Gip^~^{p) + Urip) ~ P (by virtue of the 
cancellation of the terms in G^^ and G^^~^). In this case the integral is ~ J d^p{l/p^){p) 
which is linearly divergent. This is better but still unacceptable. To ensure a finite answer 
we must choose an Ansatz which automatically enforces the right 0{p'^) and 0{p) behavior 
in G^^{p), namely, 

1 o^nBCi) 

Grip) = ^ + - — 5 — + (Ansatz starting at C»(p-^)) . (3.22) 

pZ pA 

In this case, 5G{p)/5ci < p~'^ automatically, and G^^(p) — G^j^~^{p) + ^t{p) ~ P^ ■ This is 



sufficient to ensure that the integral in Eq. (3.20) will be UV finite. The same argument 
applies to Gl and A; in each case we must build in the correct l/p^ and 1/p^ behavior of 
the propagator (or 0{p'^) and 0{p) behavior in the inverse propagator) into our Ansatz^ but 
having done so, the variations 5T/5ci will all automatically be finite (unless there are IR 
problems) . 

Applying the same reasoning to the vertices, the variation of F with respect to a generic 
coefficient di determining V gives rise to a correction of order (the phase space is explained 



in Appendix B 1 ) 



5r 

8di 



' pdp qdq kdk ^ G{p)G{q)G{k) {v - Vbare - 5^) • (3.23) 



If we allow our Ansatz to change V on the scale of its leading behavior p,k,q then this 
expression is quadratically divergent. Therefore our Ansatz must be restricted such that 
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V takes on its correct (free) asymptotic limiting behavior. Even so, in this case SV/Sdi 
and V — Vbare will be giving rise to a log divergence. Therefore we must compute 
and implement the first subleading behavior of V and only allow our Ansatz to change 

V at NNLO, 0{p~^,k~^,q~^). This will ensure finite variations of F with respect to the 
parameters di (again assuming there are no infrared issues). 

It is not necessary to determine the NNLO behavior of either self-energy or vertex cor- 
rections in order to avoid potential divergences. This is a good thing, because the (9(p~^) 
propagator correction (or 0{p^) self-energy correction) is where nonperturbative physics 
first arises. To see this, consider the one loop self-energy diagram in Eq. (3.5). Let us 
estimate the contribution when the external momentum p is large but one of the internal 
propagators is at a small momentum q ~ g"^. There is a factor of g^p^ from the vertices, 
from the hard propagator, and J d^q/q^ q g"^ from the momentum integration and soft 
propagator in the soft region. So the contribution when one propagator is soft is of order 
Hsoft ~ g"^- This contribution is nonperturbative because the behavior of the propagator at 
small momentum is. Therefore we actually cannot determine the NNLO behavior of the 
propagator at large momenta; the perturbative expansion we alluded to earlier actually fails 
at this order. Fortunately, determining this order turns out to be unnecessary to eliminate 
divergences and render our extremization problem well posed. 

Note that the elimination of divergences, both in subdiagrams and in the final variation 
of F with respect to variational Ansatz parameters, is much easier in D = 3 spacetime 
dimensions than it would be in D = 4. In that case, self-energies would be quadratically 
divergent at all loop orders and 6T/6ci would generically be quartically divergent. It is 
therefore not completely clear to us how our procedure could be extended to four dimensions 
without some changes or restrictions. We will not address this issue further at this time. 



IV. VARIATIONAL ANSATZE 

We will now start to actually solve, rather than discuss, the problem by writing out the 
variational Ansatze used for all propagators and vertices. The only dimensionful constant 
in 3D Yang-Mills theory is g. Moreover, a three-loop truncation of F only contains planar 
diagrams, and the only subleading in correction which enters when evaluating them is 
an overall factor of (A^^— 1); that is, the coupling and group theory factor for an m-loop 
bubble diagram is {N'^—l){g'^N)"^~^. Therefore, to the loop order we work, the coupling 
expansion is strictly an expansion in the t'Hooft coupling g'^N. In D = 3 dimensions the 
t'Hooft coupling has dimensions of energy, and it therefore sets the natural energy scale in 
the problem. Therefore we will factor out the overall (A^^— 1) and will scale all dimensionful 
quantities by the appropriate power of g'^N, i.e., quantities with dimension [mass]", are 
expressed in units of {g'^N)°'. For the most part, this eliminates any explicit reference to g"^ 
or in what follows. 

Variational coefficients will be generically denoted by q. This is a slight abuse of notation, 
and one should recall throughout that the q are independent for each function. For instance, 
it should not be interpreted from expressions like Gr(ci;p) and Gi,(ci;p) that Gt and 
are defined by the same set of parameters. In practice we will use rational functions (Fade 
approximants) for our variational Ansatze; we distinguish coefficients in numerators from 
those in denominators by labelling the former by a, and the latter by 6j, so that {cj} = {aj, bi} 
(or {cj} = {ttijkybijk} for vertex function coefficients). Finally, it is now always implied that 
when we refer to a correlation function, we are specifically referring to its Ansatz. We will 
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drop the q from the arguments of these functions, so G{p) imphes G{ci;p). 

Continuing on, G is first decomposed into transverse and longitudinal components 



G^-'ip) = GMT^-'ip) + GL{p)L^''ip). (4.1) 

From the arguments of the previous section, we know that whatever we write down for Gt{p) 
has to converge to Eq. (3.22) at large p (and likewise for Gl{p) and A(p)). Hence, 

Grip) = 1^^, ^tIp) = C^r^p + (Ansatz staTting at C»(/)) (4.2) 

p^ - IIt{p) 

Gl{p) = ^r^' ^l{p) = C^Pp + {Ansatz starting at £»(/)) (4.3) 

1 C^'^ 

^(P) = wTT^' ^(P)/^' = + ^t^^^i^g 0{p-^)) (4.4) 

^"^(1 — L(pj/p^j p + a; 

with 

4" = ^-^. C<"=0. C<" = 1 (4.5) 

We see that C'^pjj^i^ must be measured in units of g'^N since [C^^^^^^] = [mass] , while [Gt] = 
[mass]~^. 

The Ansatz for A differs from Gt and Gl\ namely, we assume that S(p) oc p^, and 
therefore A oc at small p. This is a condition which arises due to the structure of ghost 
vertices, which we will discuss in a little more detail when we present the ghost vertices. 
The parameter u is not treated as a variational parameter; instead its value is fixed to w = 1 
(really uj = g'^N). This choice should not be important, provided the variational Ansatz is 
flexible enough. 

We will use Fade approximants for the propagator Ansatze. The "order" of these Ansatze 
will be denoted by N^^^, which refers to the highest power of momentum appearing in the 
numerator and denominators. With this choice, Eqs. (4.2) - (4.4) read 

nT(p) = C?p + U^ip) = C?p + J^r\rr\ " . (4.6) 



i=0 "i P 



lidp) = Ci^p + UTip) = J;1f^,2^ " . (4-7) 



Z^^=0 P 



np) = c'^P + s^^(p) = ^ + J:r\,% ' ^ ■ (4.8) 

In each case we define II'^^ as the self-energy minus its one-loop perturbative (linear in 
momentum) part. In the case of E this is not the same as the part determined by the 
variational Ansatz. 

In constructing the most general gluon 3-vertex V (where it is assumed that momentum 
flows out of a vertex), six independent tensor structures require consideration. We will adopt 
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the basis used in Refs. [261 EZ] , which is 



A 


= ^mM2(Pi -P2)^3 






= ^mA*2(Pl +P2)m3 




P 


= [Pl ■ P2^/.iM2 - PlA.2P2/.i] (Pl - P2)m3 




F 


= [Pl ■ P2C/M1M2 - P1a«2P2mi] [(P2 ■ P3)PlM3 " 


- ij>i ■ P^)P2^,^] 


1 1 

■n-/Jl/^2A'3 


= ^MiM2 [(Pl ■ P3)P2,.3 - {P2 ■ P3)Pl/.3] + 3 


(P1M3P2miP3m2 - P1m2P2m3P3mi 


c 


= PIM3P2M1P3/.2 +P1M2P2m3P3mi- 





(4.9) 

Color dependence can be factored out of the vertex function, 

nr™(Pi'P2,P3) = F^^^-^W,,,,,M.P2,P3) (4.10) 

where F"'^'^ = —ifabc is the adjoint representation matrix, satisfying p^'^^pi^"- = Ca6^^, 
Ca = N for SU(A^). In a similar fashion to the propagators, our Ansatz for V is designed 
so that it can be easily made to converge to its perturbative form at large momenta. We 
will separate the various contributions to V much like we did with the propagators, 

^MlM2M3 = ■S'(^ilM2/i3 + ^ilM2/i3 + ^^1^2^*3)- (4-11) 

The bare and one-loop corrections are denoted by V^'^^ and V^^\ while V^^ denotes the 
nonperturbative correction to the vertex that has to solved for self-consistently by finding 
the stationary point of T. The bare term, V"'-'^'' is simply 

^i22/.3 = - P3)^.^g^.2^^s + {Ps - Pi)/.2^miM3 + (Pi " P2)m3^miM2 (4-12) 
which is expressed entirely in terms of cyclic permutations of the A tensor, i. e. 

V^'l,, = + cyclic perms. (4.13) 

with = 1. The one-loop correction to the vertex V^^"^ has a much more intricate tensor 
structure 

Vill,,iPuP2,P3) = A^'\P1 )P2;P3)A^iAt2/i3 

+ s«(pi,P2;P3)B^ 
+ C'^^\pi,P2;P3)C^,f,2i^, + ^^^Hpi.P2;p3)f^im2M3 

+ i7«(pi,p2,P3)H^ 1M2M3 + \Pll P2^ Ps)^ 1^1112113 

+ cyclic perms. (4-14) 

and likewise for V^^ . As we discussed above, we need the explicit forms of the 1-loop vertex 
corrections in order to eliminate logarithmic UV divergences in the variational problem. We 



present explicit results for the 1-loop vertices in Appendix A 3 and we use those results in 
the following^. Note in particular that S^^^ = 0. However it remains to write variational 
Ansatze for V^^ . Here we make no assumptions about the vanishing of the coefficients for 



Setting w*-^' = 1/4; this parameter serves an analogous purpose to uj as it appears in the ghost propagator, 



see Eq. (A35) 
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any of the tensorial structures. The functions A, C and F are symmetric in their first two 
arguments, B is antisymmetric in its first two arguments, H is fully symmetric and S is fully 
antisymmetric. We respect these symmetry properties by choosing the following Ansatze: 

Pi + P2 + Ps + ^ Y.i>j by {Plpi + PiP2)P3 

B (Pl,P2,P3) - „2 , „2 , „2 , ,.2 7^ ,{B}. i j , J i. k ^^'"^^^ 

Pi + P2 + Ps + ^ Ei>i ^Jifc (PlP2 + PlP2)P3 



^ = + + + ,{c},,.,, , (4-17) 



J Ei>i ^Iff (piP2 + PiP2)P3 

Pi +P'2+PI+^' E,>, C (^'i^'2 + PiPl)P3 
NP/ N ^ 1 E.>, 4? (^1^2 +PiP2)pt , . 

^^^'^^'^^^ p?+pi+pi+-^E.,,C(PiP^2+PlPDpt ^ ■ ^ 

i/NP. ) ^ ^ E.>,>.«g^(py2Pt + perms.) 

^^^'^^'^^^ Pt+Pi+PH-^E.,,,.C(py2P3^ + perms.) ^ " 



'7NPf„ „ „^ - - ^i>j>k%k i^^yzPxPjPzJ , . 

Pi + P2 + Ps + ^ Li>i>fe (P1P2P3 + perms.) 



which automatically have these symmetries built into them {e^yz is the permutation symbol). 
Each sum is truncated so that (i+j+k) < A^^ax- 

Notice that not all six tensors are of the same dimension; A and B have dimensions of 
[mass], C, H and S are [mass]^ and F is [mass]^. In every case the UV behavior of the 
vertex function must satisfy V'^^'^ oc 1 and V^^ oc p~^. Wc enforce the correct momentum 
scaling for the vertex functions by hand, so that the Fade Approximants are all C(l); this 
way, despite the individual dimensionalities of the Vertex functions, all of the Ansatze are 
the same "size." 

The ghost-gluon vertex V is somewhat simpler. Factoring out explicit color dependence, 

V™(pi,P2,P3) = 5i^"^"^"^V^3(pi,P2,P3) (4.21) 

with the outgoing ghost and gluon indexed by (ai,pi) and (a3,/^3,P3) respectively (pi flows 
outwards), we have 



where 



and 



dNP 



V,,3(P1,P2,P3) = A(pi,p2,P3)PlM3 +®(Pl>P2,P3)P2/.3 , (4-22) 

Hpi,P2,P3) = A(°) + A(i) + A^P (4.23) 
,,P2,P3) = B«+B^P, (4.24) 



(4.25) 
(4.26) 



1 E «iffc PiP2Pt 



NP _ ^ ijk 



pi+pi+pi + ^'Eb]j!:p\p'2Pl 

1 Pl E Q^ffpiP^Ps + (P2 - P3) E Qi?V2P3 

P?+Pi+Pi + ^' EOi^^2P3' 
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In this case A'^''^ = 1, and as above, A*^^) and B*^^^ are one-loop corrections (read directly 
from Eqs. (A36) and (A37), modulo g'^N). 

The Ansatz for B is chosen so as to guarantee that limp^_>oB = 0, linearly in pi. Let us 
briefly discuss this assumption, and our similar assumption that the ghost self-energy oc 

at small p. Both properties arise because the tree-level ghost vertex Vj?''(pi,p2,P3) oc pi^ 
the outgoing ghost momentum. In any loop diagram modifying V, with arbitrarily many 
loops, the outgoing ghost line, with momentum pi, always encounters a bare V, leading to a 
proportionality of the full diagram to pi. This proportionality is automatic in the A term; 
we are also enforcing it in the B term. One could argue that this argument assumes a strict 
diagrammatic expansion and might be violated somehow when we fully resum. However, it 
is at least self-consistent that B oc pi (note that lim^^.^ 0(^2—^3) = linearly in pi). And 
if it is, this is enough to ensure that any vertex correction, with resummed as well as bare 
vertices, is still always proportional to pi. Since the vertices are always proportional to pi, 
the self-energy must also vanish at least linearly in p; but assuming that S(p) is smooth 
at small p (which is true provided that there are no IR divergences in ghost self-energy 
corrections), S(p) must in fact vanish quadratically in p. Building these properties into our 
Ansatze improves the stability of the numerical extremization; however we have also tried to 
solve the variational problem without these assumptions (using Ansatze which allow B oc 
and S(p) oc with results which are consistent with the assumed oc pi and oc behaviors. 



A. Numerical Implementation 

Obtaining the solution in terms of the variational coefficients involves performing three 
non-trivial tasks, which together can be referred to as the numerical implementation. These 
tasks are 

• Tensor contraction and diagram generation 

• Numerical integration over a 9D phase-space (three loops) 

• Using an extremization algorithm to locate the extremum of F. 

Concerning diagram generation, the purely gluonic Mercedes-Benz is by far the most 
complicated diagram. Each propagator has 2 tensorial structures, each vertex has 14 (the 
permutations of the 6 structures described in the last subsection). An inefficient tensor 
contraction would therefore contain 2^14^ = 2458624 terms. Therefore it is important 
to perform the tensorial contractions carefully, building intermediate structures with the 
minimum number of terms. For instance, the Mercedes-Benz can be regarded as 




(4.27) 



which we will write as vertex contracted with triangle. The triangle can be represented 
in terms of a basis of 36 tensors (all three-index objects that can be constructed out of 
three momenta [two external and one loop] and the metric), and likewise, the resummed 
vertex contains 14 distinct tensors (all three-index objects that can be constructed out two 
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momenta and the metric, due to momentum conservation). Specifically, the vertex after 
contracting the tensors associated with the propagators is of form 

^l/i2Ai3 — Z00l9fllt^2PlfJ.3 + ^0105'Atl/i3PlAi2 + ^1005'm2M3P1M1 

+ Zoo2gfii^i2P2ii3 + ^0205'Mi/i3P2/^2 + ^2005'M2M3P2Ati 

+ ^112Pl/.iPlM2P2/.3 + Zl2lPl^,^P2^,2Pl^^3 + ^211P2/.iP1/.2P1M3 + ^lllPlMiPlM2PlM 

+ ^221P2/.iP2m2P1/^3 + ^212P2miP1/.2P2/.3 + Zl22Pl^,^P2^,2P2|,3 + ^222P2miP2m2P2/.3 • (4-28) 

With the vertex and triangle factored as such, standard algebraic packages can perform all 
of the remaining tensor contractions and simplifications. 

In addition to drastically simplifying diagram construction, the use of these bases allows 
for an economical use of floating point operations. The triangle is contained in all gluonic 
derivatives of F; furthermore; the 36 triangle Z-coefficients are by far the largest polynomials 
contained within the problem. With the triangle expressed as such, each of these 36 functions 
need only to be computed a single time at each point in a 6D space of integration variables. 
This is important because high-dimensional numerical integration will require of order 10^ 
evaluations of each diagram per step in the extremization procedure for F. 

With the Lorentz algebra in hand, we turn to the problem of multidimensional numerical 
integration. The first step is to choose a convenient basis for integration. Our choice 
is described in Appendix |B] Performing the global (Eulerian) angular integrations, two 
loop diagrams require a 3D integration and three loop integrals require a 6D integration. 
These happen to be the same as the number of propagators in the bubble diagrams built 
with 3-point vertices. And it is possible, and convenient, to choose integration variables 
which are precisely the magnitudes of the momenta on each propagator. (This is a special 
feature of phase space integration in 3 dimensions, discussed in the appendix.) The most 
numerically challenging integration is again the Mercedes-Benz topology. In the notation 
of Appendix [B| the three finite-range (angular) integrations, over k', q', and /, as well as 
the p integration, were performed using Gaussian quadratures, while the (infinite) k and q 
momentum axes were rescaled into the unit interval and sampled using an array of points 
constructed with a quasi-random hopping Halton series.^ Our quadratures procedure is 
symmetric over intervals, which ensures that certain cancellations on angular integration are 
preserved, and it avoids edgepoint evaluations. Also, since neither algorithm uses random 
or pseudorandom numbers or dynamic mesh refinement, each integration evaluation is over 
exactly the same distribution of phase-space points. This ensures that the effective action 
F is not "noisy" in the sense that it does not fluctuate between evaluations with the same 
or almost the same choices of propagator and vertex functions, a feature which is essential 
for conjugate-gradient and other differential extremum seeking algorithms. 

To test the stability of our algorithm against changes in the number of integration points, 
we computed the 2-loop self-energy for particular values of full propagators and vertices with 
varying numbers of integration points. As illustrated in Fig. |3} the results converge when 
sufficiently many integration points are used. 

Now that we have explained how numerical integrations can be performed, we turn to 
the problem of extremizing the effective action F. One challenge is that, because of gauge 



^ More specifically, we first write k = px/{l — x) with x G [0, 1). Tlien we write x — — 2y'^ with y e [0, 1) 
chosen with uniform weight. The former transformation ensures that fc, q,p are of comparable magnitude; 
the latter transformation increases the sampling at the top and bottom relative to the middle of the range. 
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FIG. 3: Non-perturbative correction to 11^. viz. Eq. (4.29), as a function of the total number of 



integration points along the fc, q, k' , q' and / integration axes (labelled in the bottom right corner). 



fixing and ghosts, the extremum is actually a saddle, rather than a maximum or minimum. 



This can be easily seen from the first line of Eq. (2.6), where the gluon propagator G and 



the ghost propagator A enter with opposite sign. With the sign conventions chosen there, 
the extremum of the 1-loop action with respect to G is a maximum; with respect to A 
it is a minimum. This rules out any straightforward application of the conjugate gradient 
algorithm. The Newton- Raphson algorithm can find general extrema, but it is inefficient and 
tends to converge well only in rather small basins of attraction. So some hybrid approach 
is needed. 

Fortunately, despite it being not at all a priori obvious, one observes for the most part 
that each individual function that makes up the extremization problem has relatively little 
effect on the others. This opens up the possibility of iteratively extremizing each constituent 
function (G, A, V and V). Our procedure was to start with vertices set to their 1-loop values 
and propagators set to some naive initial guess. Then we perform a conjugate gradient 
extremization with respect to G, then a conjugate gradient extremization with respect to 
A (with opposite sign on the gradient).^ Iteratively extremizing G and A solves the 3-loop 
2PI problem. Then we use gradient descent to extremize F with respect to the three-gluon 
vertex functions. Finally, the ghost vertex corrections are improved using the Newton- 
Raphson method. Then the procedure is iterated (propagators and gluon vertices, then 
ghost vertices) until convergence is achieved. 

The convergence of the algorithm, when applied to the 3-loop 2PI problem in Landau 
gauge, is depicted in Fig. El The figure compares the Ansatz value for the self-energies, 11^^ 



In practice we also accelerate this procedure as follows. We evaluate the self-energy diagrams n(p) at a 



sample of p values holding G fixed. Then we conjugate-gradient extremize G in Eq. (3.201 but treating 



the self-energy 11 as fixed. We insert the new value of G into the evaluation of the self-energy and iterate. 
This minimizes the number of evaluations of n(p), the most numerically expensive part of the procedure, 
needed to converge to the extremum. But the extremum obtained by this procedure is the one which 



satisfies Eq. (3.20), as desired 
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FIG. 4: Evolution of the self-energy under the gradient descent algorithm. Here is the 

nonperturbative self-energy according to the Ansatz, while 11^^ is the value as determined by 
evaluating the self-energy diagrams. As the algorithm is iterated, the Ansatz approaches a correct 
reproduction of the self-energy. 



and E^^ as defined in Eq. (4.6), Eq. (4.8), to the values directly evaluated by summing 
the self-energy diagrams, and S^^; in each case we have removed the 1-loop linear- in-p 
contribution. That is. 



^T^^(n2 + ng)/(D 



(4.29) 



and similarly for S^^. The extremization procedure for the effective action with respect to 



one of the variational coefficients in Gt, Eq. (3.20), then corresponds to 



6r {D-1){N'^-1) f 6Gt{p) 



6cj 



(2n) 



6c4 



(-n?P(p) + n|°(p)) 



(4.30) 



Fig. |4] shows two things. First, even though the initial guess for the self-energy falls 
quite far off the actual value, after relatively few iterations the fitted and true values of 
the self-energy become similar, and the eventual convergence is excellent. Second, the value 
of the self-energy U^P actually depends quite weakly on precise form of 11^^ (p). That is 
why our procedure of varying F with respect to individual functions (rather than trying 
to do everything at once) works so effectively. The presence of vertices definitely makes 
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matters more complicated; however, it is also observed that Il^{p) is fairly insensitive to 
their inclusion. 



Eq. (4.30 ) should be interpreted as an Euler-Lagrange type of equation for 11, and a vertex 
analogue can be defined as follows. For instance when Cj belongs to the gluon if-function, 
the variation of F takes on the following form 



5cj 



I - / 1^(0^^^ (-H. ."^(p,,.) .H. .^o,,,,.)) , (4.31, 

with 



where the overall minus sign comes from the ordering of color indices. H ■ V^^{p,k,x) is 



defined in a similar manner, except that the term in brackets in Eq. (4.32) contains all of 
the higher loop terms in the vertex Schwinger-Dyson equation, Eq. (2.9). The values of 
H • V^^{p, k,x) and H ■ V^^{p, k,x) along the curve defined by = p/4 and cos 9pk = 1/4 
are plotted in Fig. |5| Fig. |6] contains plots of a similarly defined set of quantities related to 
the ghost vertex. 

The Ansatze in Figs. [4]- [g] correspond to A^j^ax ~ 3 and A^j^ax ~ 3, from which we observe 
that the solution is well described by third order Fades. However, the size of the Ansatze can 
have a major effect on the outcome of this technique. With too few coefficients, the numerics 
are much simpler, but one is not able to obtain the correct final answer. However, as the 
number of coefficients increases, the problem becomes very numerically difficult (considering 
the relative ease with which poles may form in the denominators of the Fade approximants), 
and furthermore, the data will eventually be over-fitted. We find that at A^^ax = 3 and 
-^max = 3 the problem is still not overly difficult to solve, despite there being a total of 174 
coefficients in Landau Gauge (with an additional 33 in other gauges). At the same time we 
are not over-fitting. This choice is further motivated by the general shape of functions we 
are attempting to converge to; the presence of additional "wiggles" would necessitate larger 
Ansatze. 



V. RESULTS 



The extremization procedure was carried out for several choices of the gauge parameter 
^, namely 0.0, 0.5, 1.0 and 2.0. The resulting self-energies are shown in Fig. [7], and the gluon 
three- vertex functions A through S as well as A^^ and B^^ from the ghost-gluon vertex are 



plotted Figs. [12] and [13] for Landau and Feynman gauge respectively. The Landau gauge 
variational coefficients are stated in Table |TTj It is observed that the solutions in ^ = 0.5 and 
^ = 2.0 are qualitatively similar to what has been depicted so far. 

In Landau gauge, Gl is zero, so it is not included in the variation; hence Gl and H^ are 
depicted as zero in the plots. Furthermore, when the tensorial structure B^^^j^g is contracted 
against transverse propagators on all three legs, the result is zero; therefore the coefficient 
B vanishes exactly in Landau gauge, though not in gauges with nonzero ^. The function 5" 
turns out to vanish in all gauges. 

The effect of including the vertices and allowing them to vary is shown in Fig. [8] where 
we see a comparison between the 2PI and 3PI solutions. The inclusion of the vertices only 
has a slight effect the resulting propagators. 
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FIG. 5: Vertex analog of Fig. [4] in Landau gauge, illustrating in this case, convergence of the H- 
and cos 6pk 



function (see Eq. (4.31)) and the gluon vertex as a whole. The values reside on the curve k 
1/4. The figure on the left corresponds to H^^ 



0, whereas the figure on the right 
corresponds to H^^ which extremizes F. Note that the basis, Eq. (4.9), used for the vertex is not 
orthogonal, so H^^ = does not imply that H • ^ 0, as illustrated here. 
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FIG. 6: Convergence of the ghost-gluon vertex functions in Landau gauge; similar to Fig. [5| the 
figure on the left corresponds to an initial guess of V^^ = 0, whereas that on the right corresponds 
to the "solution" for which extremizes F. 
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FIG. 7: Self-energies. 
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p p p 

FIG. 8: Comparison between the 2PI (bare vertices only) and 3PI (vertices included) solutions. 
The correction to Gt that we obtain when including the vertices is indeed small (in Landau gauge). 

The dependence of G, V , etc. on the choice of the gauge fixing parameter ^ does not 
by itself indicate a breakdown or hmitation on the 3PI approach. The relevant question is, 
how dependent are gauge independent quantities on and how closely do such quantities 
correspond to the nonperturbative values determined, for instance, using lattice techniques? 
Any ^ dependence in gauge invariant quantities would be an ambiguity, and any error in 
their value in comparison to lattice determinations would be a failure, of the 3PI technique. 
Such comparisons are essential, but they are beyond the scope of the present manuscript. 



VI. DISCUSSION 



A. Comparison with other approaches 

The majority of the literature on this subject is centered around 4D Yang-Mills theory; 
however, lattice studies (described shortly) have shown that Green's functions in 3D and 4D 
exhibit similar qualitative behavior. Nevertheless, we will try our best to directly compare 
our results with those obtained in 3D, to the extent that they exist. 

The gluon propagator is an interesting quantity, despite not being directly related to 
any physical observable. Gt as depicted in Fig. [8] violates reflection positivity, that is, 
it does not have a Kallen-Lehmann representation in terms of a positive spectral density. 
Hence, in (3+1) dimensions (or in our case (2+1) dimensions) Gt can not describe the 
correlations of physical particles. This violation of reflection positivity is allowed despite 
its apparent contradiction with the Osterwalder-Schrader axioms; after all, we are dealing 
with a confining theory, so there is no one-to-one correspondence between fields and physical 
particles. This is further discussed in greater detail in [28], but the main point is that this 
behavior signals confinement. 

The only propagating degrees of freedom that we can (in principle) observe are color 
singlet bound states, gluehalls for instance. Hence, gluonic two and three point functions 
are not "physical," and, in general the results we have presented are ^ dependent, which 
is not necessarily a bad thing. Indeed our intention is to use these results to compute 
gauge-invariant observables in some later publication. 

However, in the mean time, we are stuck with analyzing {AA) and (cc) as well as the 
vertices. Though arguably our best insight into IR QCD comes from the lattice, there have 
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been many notable first-principles based speculations about the specific IR form of these 
functions. In general, the main point of contention is the exact value of Gt{0)- The most 
popular schools of thought can be summarized as follows... 

1. The Gribov-Zwanziger Confinement Hypothesis 

In his study on gauge-fixing and gauge copies in Yang-Mills theory [22], Gribov proposed 
that the IR behavior of the gluon and ghost propagators is 

Gt(p)~^^, A(p)|,.^o~^ (6.1) 

which has the generic feature that Gt vanishes at zero momentum, and moreover, A experi- 
ences IR enhancement, which can possibly be interpreted as signaling linear confinement 
(in 4D, or course). This form of Gt and A was later advocated by Zwanziger, primarily 
because it vanishes at p = 0, which is in accordance with his theorem that Gt(0) (in Landau 
gauge) must vanish on any finitely spaced lattice in the infinite volume limit |30j . 

This proposal should be regarded as being fairly dated, and it is not in agreement with 
any of the more recent lattice data. It also does not accord with the behavior we determine 
by solving the 3PI problem. 

2. Schwinger-Dyson Equations 

These arguments [31] are based on obtaining solutions to a truncated set of Schwinger- 
Dyson equations for the gluon and ghost propagators, and in a sense, are very reminiscent 
of what we are doing here. If we assume power behavior of the gluon and ghost propagators 
in the infrared, 

p'Grip) ~ (p'r-, p'A{p) ~ (p2)-«- , (6.2) 
then Ref. [22] claims that kq, i^a must satisfy 

/€G = 2/€A + (4-D)/2, (6.3) 

and specifically in 3 dimensions = 0.2952, implying that the gauge field propagator goes 
to zero and the ghost propagator diverges more strongly than l/p"^. However this result 
assumes that the loop integral giving rise to a self-energy at momentum p is dominated by 
momenta of order p, whereas we find for small p that it is instead dominated by momenta of 
order g'^N. Therefore it is not clear to us that this result of Ref. [22] is robust, see also [H3] . 
It is also contradicted my more recent studies [SUES], which give results (in 4 dimensions) 
showing Gt{p) going over to a constant, and A(p) oc in the infrared. These studies 
are in at least qualitative agreement with lattice investigations. However, since in general 
Schwinger-Dyson based approaches are reliant on many simplifying approximations, they 
have yet to produce any quantitative agreement. 

3. Observations From the Lattice 

There is a wealth of lattice data related to this subject, and fortunately, different sources 
are generally in agreement. Simulations have been performed on very large lattices {V = 96^ 
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FIG. 9: A comparison between the results of our calculation {Gt{p) in Landau gauge) and a 
reproduction of the plot of aDijp) (appropriately rescaled) in [|iDJ. A free propagator is shown 
for reference. 

V = V = 128^ and V = 320^ [3H]), from which one observes quahtative 

agreement between the results for 3D and 4D (hence we will intermittently compare 3D and 
4D data, but never 2D). The generic finding is that p'^^{p)\p2^o and Gt{0) are finite and 
nonzero. 

However, Gt{0) is often seen to scale inversely with volume so it remains an open question 
as to whether the Zwanziger hypothesis is observable, and it is not known at what volume 
one should expect to see this effect. The results in Landau gauge currently depict a l/V^ 
scaling for Gt but it is generally not observed that Gt{0) — )■ as \^ — )■ oo. 

All of the works cited above specifically employ the lattice implementation of Landau 
gauge. In fact, it is only fairly recently that preliminary 3D and 4D results in Feynman 
gauge have been made available 

In Fig. |9] our data is compared directly to the results in [IQ] . Their calculation was 
performed for SU(2) on an 80^ 3D lattice with /3 = A/ag'^ = 6.0. To facilitate the comparison 
we have recast their results so momentum is scaled by g'^N. In the infrared their results 
are qualitatively similar to ours but differ quantitatively at the factor-of-2 level. This might 
indicate a limitation of the large- iV expansion, or it might simply indicate a failure of the 
3PI method. 

B. Slavnov- Taylor Identities 

Planar diagrams on their own form a gauge invariant subset of the full loop expansion 
[IT] . One may hope that in resumming a "dominant" or "important" set of planar diagrams 
(which is what hope to be doing here) gauge invariance is approximately conserved. This 
can be measured seeing to what extent our resulting two and three-point functions violate 
the Ward-Slavnov- Taylor (WST) identities. For the gluon propagator, we have 




(6.4) 



with deviations from this identity shown in Fig. 10 
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FIG. 10: Propagator Ward Identity 




FIG. 11: Vertex Ward Identity 
In its most general form, the WST identity for the gluon three-vertex is 

^ {pl9a^^2-P2aP2,,2)^'^,{P^,PuP3), (6-5) 



JiP2) 



where F and J are defined in this context as J{p) = p'^Gt{p) and F{p) = p'^A{p). V with 
two Lorentz indices is given by 

V°3(pi,P2,P3) = 9'^:,a{P3,P2,Pl) - P3P2ti3b{P3,P2,Pl) + PiP3tiAP3,P2,Pl) 

+ PtPit^zd{P3, P2,Pi)+ PiPiiise{p3, P2,Pi) (6.6) 

which is related to the usual ghost-gluon vertex via V^3(pi,p25P3) = Pia^'^^iPi, P2, P3) (fol- 
lowing once again with the decomposition in [26]). From this identity one obtains 

pfprpr^^,M2M3(pi,P2,P3) = o. (6.7) 
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With p'^ = p^/p, the deviation from the vertex ward identity is show in Fig. 11 for Landau 
gauge. As previously, the vertex is a function of 3 variables, so to make a 2+1 dimensional 
plot we have fixed an angular variable to cos6'p^p2 = 1/4. 



VII. CONCLUSIONS 

We have successfully found the propagators and vertices which extremize the 3-loop, 3- 
particle- irreducible action of QCD in 3 dimensions (the all-threes problem). We did so by 
writing a nonlinear variational Ansatz for three propagators (ghosts and the transverse and 
longitudinal components of gluons) and for eight vertices (two tensor structures for ghost- 
gluon vertices and six tensor structures for three-gluon vertices). To avoid divergences it 
was necessary to add and subtract terms to 1 and 2 loop self-energies; the added terms 
are computed in MS, the subtraction renders the remaining numerical integrals finite. It 
was also necessary to compute the first loop corrections to 2-point and 3-point functions at 
large momentum explicitly and to incorporate these corrections into our Ansatze for those 
functions. 

The most urgent task is to test the resulting resummation against exact nonperturbative 
results in 3D QCD by comparing the values of gauge invariant questions. For instance, one 
should be able to evaluate the (F^(x)F^(O)) correlator, whose Fourier transform gives the 
lowest glueball mass. It might also be possible to evaluate the correlator of two field strengths 
connected by a Wilson line, which is of interest in evaluating the Debye screening mass in 
full QCD |l2]. Slightly extending our treatment to include a fundamental representation 
scalar, it should be possible to explore the 3D SU(2)+Higgs phase diagram, which can also 
be found nonperturbatively on the lattice [1]. 

Unfortunately it is not possible to compute the pressure of 3D Yang-Mills theory at the 
3PI level, because the nontrivial contributions to the pressure arise at 4 loops. Evaluating 
the pressure would require a solution to the 4-loop 3PI or 4PI problem. Extending our 
approach to the 4-loop 4PI treatment would not raise any new conceptual issues, since all 
potential UV divergences in the extremization procedure are already encountered at the level 
of the 3-loop 3PI problem. It would be interesting to do so because the nonperturbative 
contribution to the pressure of 3D Yang-Mills theory is needed to compute the term in 
the pressure for full QCD [161 US]- However the extension to 4 loops and 4PI would be 
prohibitively difficult because the diagram generation and loop integration would become 
even more cumbersome and the number of possible tensor structures for the 4-point function 
is large. 

Having studied the 3D theory, it appears to us that the extension to 4 or 3+1 dimensions 
will be extremely difficult. The problem is that one is simultaneously solving nonperturbative 
infrared physics and (perturbative) ultraviolet physics. The effective action is extremely 
sensitive to the ultraviolet form of the propagators and vertices; a procedure along the lines 
of what we have done here encounters quadratic UV divergences at every loop order when 
evaluating self-energies, and quartic divergences when varying the propagators in a way 
which changes their UV behavior - at every loop order. It would be much harder to "cover 
up" gauge non-invariance in 4D because divergently large gauge boson masses would arise 
at every loop order, whereas we only encountered them at one loop (at two loops there were 
logs but they all cancel). Thus it is not clear what additional techniques would have to 
be developed to successfully extend our procedure to 4 dimensions. We will leave this for 
future investigation after we know whether the 3-loop treatment is successful in describing 
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nonperturbative physics. 
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Appendix A: Computations involving Bare Diagrams 

Since we arc working in three Euclidean dimensions, it is useful to now state the bare 
Feynman rules (with vjf^Zt' = F-^---W^.%,,, and vgfe^f^ = i^"^"^"^vfi,^3) 

+ {P3 ) (Al) 

^Vf (pi,p2,P3) = gpi^s (A2) 

+ F-^^F'^'ig^^g^r - g^pg^r)) (A3) 

noting especially the overall sign of the four vertex. F"'"^ are the generators of the adjoint 
representation, defined in terms of the usual SU(iV) structure factors by F"'^'^ = — i/"*"^. 
piahpjba _ qj^(\ for SU(A^), Ca = N. The bare propagators are simply given by 

G'!S(P) = ^(T^^ + eV) (A4) 
A(°)(p) = 1 (A5) 

with the transverse and longitudinal polarization tensors as defined earlier (color indices are 
suppressed) . Finally the loop integral in D = 3 + 2e dimensions is defined as 



1. Gluon Self-Energy 

The gluon self-energy contains diagrams that superficially diverge linearly (one-loop) and 
logarithmically (two-loops); however, when computed in dimensional regularization, one 
finds that the sum of all one and two-loop contributions to the self-energy is UV finite. 

At one loop, there are three diagrams, only two of which are nonzero in dimensional 
regularization. 
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Keeping terms up to 0{e) for convenience, the result is found to be 



„2 AT ^l + 2t 



p 



2e 



2e 



]-fB(l,e)rp 



1X1/ 



(A9) 



UV divergences only start to arise at the two-loop level, and on a diagram by diagram basis, 
these are all proportional to g^i,. The two-loop self-energy also contains diagrams that are 
potentially IR divergent, and these will be handled separately. Defining 



^-(D — 



(AlO) 



where it should be noted that the "figure eight" diagram (proportional to 1/4) is strictly 
finite, but is included for completeness. An actual computation yields 



TTB(2,e)UV ^ 9 P 
V ^2 ,, 



(e + 2)(e^ + 2e + i) 

768e 



^ 8(7e^ + 75^^ + 22ie + 233) - 18C(2)(e^ + 3){e + 2^^ + 17) ^ 



12288 



7f + 32(- + 79(+ 12 ^ 
768 



(All) 



The IR regulation is achieved via the introduction of a fictitious mass in the denominators 
of both of the divergent diagrams. In D = 3+2e, the IR regulated integrals have the following 
form 



V. 



(0) 



(A12) 



"0^ = 



TTB(l,e) o"/5 
,2 / y(0) y{0) 11 ^ ' ^ y g-+m 



_]_')■>■) 2 iJ 



(p+9)^(p+g)'' 



(A13) 
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Including the IR regulator w? in exactly this manner was found to yield relatively simple 
and compact final expressions. Adding together these two diagrams, 




(A14) 



we obtain 



„4 at2 At 
TTB(2,e)IR Reg ^ 9 -'^ P_ 



768e 



+ ^ im + 2)(e' + 2e + 11) log ^ + ^(9^' + + 93^ - 33) 



2304 



m 



+ + 48^ + 54^ + 660) + (IS^'^ + 148^ + 299^ + 830) 



3(^2^2^ + 11) 



(3^ - l)m^ + (11 - 30 (mV + ^ P ) + (2C - 2)p' 

p4(p2 _|_ j^2'J 



+ 



tanh-^ ./^ 
V p + 



p3(j)2 _|_ 777,2^/2 



((4^ + 8)m^ + (9^ + 43)p2m^ + (4^ + 34)p'^m^) 



1 f 16?)^ TTT-^ 

+ 7^ ( 3(e + 2)(f + 2e + 11) log - ^(24^3 + 96f + 360^ + 528) 



2304 

+ (21^=^ + lOOf + 245^ + 170) 
tanh~^ 

+ 6(e' + 2^ + 11) 



p 



pS^p2 _|_ ^2^1/2 



((4^ + 8)m^ 



(A15) 



As expected all of the UV divergences cancel between diagrams, and hence we can safely 
take the e ^ limit and be left with something finite. 



hm (n^(2,^)UV + Me)m RegA ^ 



(A16) 



To be able to ensure that G reproduces the correct subleading 0{p ^) behaviour, we simply 
need to know the one-loop gluon self-energy. 11^*^^) as referenced earlier in this paper is given 

by 

^2nB(i)^n'^(i,o)_ 
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2. Ghost Self-Energy 



Knowledge of the one-loop ghost self-energy is necessary to guarantee the correct 0{p ^) 
UV behavior of A. There is only a single diagram, 

= __(/^__ (A18) 



which in D = 3 is independent of ^ and equals 

S^(i'°)(p)=p<72S^W=p^2^. (A19) 

All of the bare diagrams that contribute at Olg"^) are UV finite, so there is no need to 
carry out this calculation to the next order. 



3. Three-Gluon and Ghost-Gluon Vertices 

In this section we will present our results for the one-loop corrections to the three-gluon 
and ghost-gluon vertices valid for arbitrary covariant gauge in 3D. The generalization to 
arbitrary D is available in the literature |27] . 

When working in covariant gauge, one encounters scalar integrals of the form 



~ / TTTZvnT, ,„nrvr, .„.FSr „n -v (A20) 



(Pq 1 

where a, P and 7 can take on values between —3 and 2 over the course of the calculation. 
In three dimensions, all of the triangle integrals with integer a, (3 and 7 are finite in DR. 
Moreover, they can all be expressed in terms of 

Ctn = ^ (A22) 
Cf„. = ^ (A23) 
C?,„ = ^ (A24) 

with P3 = (pi + ^2)^ and Pi = The relations between triangle integrals with different 

a, P and 7 can be obtained from the generic expression for C°^^ which is known in terms 
of Appel's hypergeometric function 



00 00 
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making use of the Pockhammer symbol (a)j = r(a + i)/r(a). In any (Euclidean) dimension, 
omitting the MS scale, the expression reads 



C\ 



1 



«/37 (47r)D/2r(7)r(/3)r(a)r(D - 7 - /3 - a 
(p2)D/2-7-/3-ar(^)r (t + + « " D/2) T (D/2 " 7 - /?) T (D/2 - 7 - a 

(2 2 
7, 7 + /3 + a - D/2; 7 + /3 - D/2 + 1, 7 + a - D/2 + 1 ^, ^ 

+ {vlT'^-'-'^ivlr^m^ (D/2 - «) r (D/2 - 7-/3)r(7 + a- D/2) 

(2 2 \ 
/3, D/2 - a; 7 + /3 - D/2 + 1, D/2 - 7 - a + 1 ^, ^ ) 

+ (p?)°/2-^-/'(p^)-r(a)r (D/2 - /3) r (D/2 - 7-«)r(7 + /3 - D/2) 
X F4 ( a, D/2 - /3; D/2 - 7 - /3 + 1, 7 + a - D/2 + 1 



Pi Pf 

2 ' 2 

Ps Pi 



+ r(D - 7 - /3 - a)(j9?)°^'"^"^(P2)°/'"^""(P3)""°^' 
X r (D/2 - 7) r (7 + /3 - D/2) r (7 + a - D/2) 

X F4 (d - 7 - /3 - a, D/2 - 7; D/2 - 7 - /3 + 1, D/2 - 7 - a + 1 







2 


Pi) . 



(A26) 

At the one loop level, the three-gluon and ghost-gluon vertices have the following form 

^?<tMr(Pi'P2,P3) = ^7F'^^'^^'^«[A«(Pi,P2;P3)A^,^,^3 +i?(^)(pi,P2;P3)B^,^,^3 

+ C«(Pi,P2;P3)C^iM2/.3 + ^^'Hpi,P2;p3)F, 

+ /7(^)(pi,P2,P3)H/,,^2^3 + ^(^)(pi,p2,P3)S^ 

+ cyclic perms. ] (A27) 

9^'^il^'^'^'{Pi,P2,Pz) = gF'^^'"'''-^' [a^'\pi,P2,Pz)pi^, + &^'Hpi,P2,P3)P2m3] (A28) 

with the tensors A through S as defined in Section [TV] The vertex functions are as follows: 



VlA'2M3 



A^^\pi,P2]P^) 



16piP2(Pl +P2+ P3) [4(pi - P2)^(Pl + P2) 



1024pfp|p3(pi + P2 + P3)^ 
+ {5pj + QpiP2 + 5^2)^3 + 3(pi + ^2)^3 + 6P3] - 4[(p? - P2)^(Pi + 2p?P2 
+ 4piP2 + 2piP2 + P2) - 2piP2(pi + P2){pt - 4:pIp2 - 2p\pl - Apipl + p\)p^ 

- {p\ - uptpl - iqpIpI - HpIpI + P2)P3 + ^pIpHpi + ^2)^3 

- {Pl - P2f{pl + 4piP2 + ^2)^3 + 2piP2(Pl + P2)P3 + (Pl + ^2)^3] (1 - 

+ (Pl + P2 + Pzf [{Pl - P2f{,p\ + 2piP2 + 2piP2 + P2) - 2(pi + P2)^(P? - 3]9iP2 + ^2)^8 



+ 2(p1 - pIpI + P^)P^ - 2{p\ + p^)p3 + (p2 ^ ^2)^4] _ 



(A29) 
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1QpIpIpI[2{pi + P2f 



1024pfp|p|(pi +P2+P3y 
- 9{Pi+P2Yp3 - 20{pi+P2)pI - QPs] - 4[2p2(p, - ^2)^2(^1 + 



+ ¥ipi{pi - pifp3 - {pi + P2){p\ + 2piP2 - ^p\pi - ^pIpI - "ipivl 



2pipl + pDpI 



+ 2piP2{pi + P2) {Pi - QpiP2 + P2)P3 + (Pl + P2){Pl - 26piP2 + P2)P3 " ^'^PlP2P3 
+ {Pl + P2fpl - 2piP2pl - {pi + P2)pI] (C - 1) 

+ {P1 +P2+ P3f [2pi{pi - P2fpl{pi + P2) - {Pl + P2){pt + Pt)pl 

+ 2(pi - p2f{pl + 3piP2 + pDpI - 2{pi + P2){pI + pDpI 



+ 2(p? +PiP2+pI)pI - {Pi +P2)pI] (1 - 0' 



(A30) 



C^'\pi,P2;P3) = 



48piP2(Pi +P2) P3 



512pfpl{pi +P2)P3{P1 +P2 +P3y 
+ [(Pl -P2f{pi +P2f - 2pIp2{Pi +P2)^P3 - 2pipl(pi + P2)^P3 - ptiPl + P2)p 
- i+Pl + P2)pipl - {Pl + P2)V3 + 2PIP2PI + {Pl + P2)pI] - 1)(3 + 



+ 2pi{p, + P2)pipiisi + 5[2 + + ^Piplpli29 + [6 + m 



(A31) 



F^'\puP2;P3) 



[- (Pl -P2f{Pl +P2f 



^12pIpI{pi + ^2)^3(^1 +P2+P3f 

- ^{Pl - P2)\P1 + P2?P3] (e - 1)'(3 + + [iPl + P2)pI{-3 + 0(3 + 

+ (3 + e)(-20p3p2(pi + P2)'P3 - 20pip^(p, + ^2)^3 - pUpi + P2)pI{9 + 

- ptiPl +P2)pI{9 + +PiP3(ll + + 4p?P2P^(ll + + 4pipiP3(ll + 
+ P2P3(ll + -pt(Pi+P2)'pI(7 + 30 -p^(pi+P2)'pl(7 + 30)](e- 1) 
+ 2p?p^(pi +P2)V3{ - 113 + 3e(-l + [-5 + 

+ 2pIpI{p, + P2)P^( - 303 + e(-25 + [23 + 

+ 2p\pIpI ( - 215 + e(33 + [35 + 3^10) + [Kpi + P2)'pl(e - 1) 
+ 12(pi + p2)(p? + P2)P3 - pI + 14piP2 + ^pI 

-3(pi+p2re]l(l-0(-3-0l (A32) 
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+ 6 [pIp2P3 - pIpIp3 - pIpIps - p\p2pI 



1024pfp|p|(pi +P2+P3Y 

- pIP2pV\ + 3pi(p2 + Ps) + 3piP2(P2 - P3f{P2 + Ps)^ + ^PipI{P2 - P^fp^{P2 + Ps)^ 
+ QPIP2{P2 - P3fpl{p2 + Pzf + 3pi(p2 -P3)^P3(P2 +P3)^ +P2(P2 -^3)^(^2+^3)^ 
+ 2p^(P2 - P3fp3{p2 + P3)^ + 2P2{P2 " P3)^P3(P2 + ^3)^ 

+ iP2 - P3?pI{P2 + Pzf] - 1)'(3 + + [[2p?P^P3(P2 + Ps) 

+ '^pIp2pI{P2 + Ps)] (13 -0 + 2 [pIpI + + pIpI{P2 + Ps) + P^Psfe + Ps)] (7 - 
+ [2pIpI - 4p^p^ + 2pIpI - Ap\pl - Aplplip^ + p^) + 2p2p6(p2 + 

- ^p\pI{P2 + P3) + 2piP3(p2 + Ps) + 6piP2(P2 - P3)^P3(P2 + ^3)^ 
+ 2pl{p2 - P3fpl{P2 + P3f] (5 + + 2pW3(P2 + P3)(ll + 

- [2ptp3p3(p^ + + 2ptp2pliP2 + P3)] (23 + 0] (e - 1)(3 + 

+ MpIpI{P2 + P3) (78 - + 7e=^) - 2 ^ p3^2^4] ( _ 225 + ^(53 + [25 - 13^10) 

- 4p?p^p3( _ 217 + 3^(17 + [7 _ 5^]^)) + Aplplpl{ - 18 + ^(13 + [20 + ^]0) 
+ WiP2pI{P2 + Pz) + 2pIpIpI{p2 + P3)\ (3 + e(-3 + [29 + 3^]^) 



+ ^p\pIpI{P2 + P3) (111 + ^(-25 + [-27 + 5^10) 

S^'\PUP2^P3)=^ 



(A33) 



(A34) 



and, for instance, the regulator w'-^-' mentioned in Section IV is included in the A function 
by making the transformation 



1 



1 



PfplPsiPl +P2+P3y ' {Pl + U;(l))3(p2 + W(l))3(p3 + + P2 + P3 + UJ^^^y 

and likewise for B through H. The one-loop ghost vertices are 



(A35) 



0'''^\Pl,P2,P3) 



lQpl[-pl+pli-p2+P3) 



512pip2pl{pi +P2+P3) 

- {P2 - P3f{P2 + P3) + Pi (3^2 + 2P2P3 + P3)] + 4(P1 - P2 - P3) [{p\ " plf 

+ 2pi{pi - P2)(P1 + P2)P3 + (5pi - P2)(Pl + P2)P3 + 2(Pl + P2)P3 ^ 2p^] (1 - Cl 

- (Pl - P2 - Ps) [(Pl - P2)^(Pl + P2) + (Pl - P2)^P3 + (Pl + P2)P3 - 3P3](1 - 



(A36) 



b^^\pi,P2,P3) 



g'N 



32pipl[{pi -P2)P2 +P3] 



512piP2P3(Pl +P2 +P3) 
- 4(Pi - P2 + P3) P? + (P2 - P3) [(P2 + P3)^ - Pl(P2 + 3p3) - pI] (1 - 

+ (P? - P2 + P3) [iPl - P2f{Pl + P2) + (Pl - P2) V3 + (Pl + P2)P3 - 3P3] (1 - 0' 



(A37) 
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Appendix B: Phase space integrations 

In performing numerical integrals over vacuum diagrams encountered in this paper we 
need efficient parametrizations of the phase space integrals. At two-loops the most interest- 
ing diagram is the Setting Sun, at three-loops it is the Mercedes diagram. All other diagrams 
can be solved by reducing these two (as we will describe), so we will concentrate on them. 

In D dimensions an n-loop diagram involves nD real integrations. However the symmetry 
group 0(D) helps reduce this because certain angular integrations are trivial. Namely, there 
are D(D — 1) /2 global angular integrations. Selecting n D-vectors reduces 0(D) to 0(D — n) 
(for n < D — 2) or reduces it completely (for n > D — 1). Therefore, for n < D — 2, 
D(D — l)/2 — (D — n)(D — n — l)/2 = nD — n{n + l)/2 of the integrals are global angular 
integrals which can be performed immediately since none of the invariants depend on them. 
This leaves n{n — l)/2 nontrivial integrations, for n < D — 2. For n > D — 1 there are 
nD — D(D — l)/2 nontrivial integrations. 

In an n-loop connected vacuum diagram built entirely with 3-point vertices there are 
3n — 3 propagators. For D = 3 and n > 2 this happens to equal the number of integration 
variables. Therefore, in D = 3 dimensions, in diagrams composed using 3-point vertices and 
where each propagator has a distinct momentum (which is the case for 2PI or 3PI diagrams), 
it should be possible to arrange for the integration variables to be precisely the magnitudes 
of the momenta on all propagators. This is a very convenient choice, provided that all dot 
products of propagator momenta have simple enough expressions. 



1. Two loops: Setting Sun 

We apply these ideas first to the Setting Sun diagram, that is, two vertices connected by 
three lines: 

k 




q 

The "natural" integration variables are 

/ w = r "^''^ r 1^"°"'" 

where we have performed the trivial integral over the Eulerian angles, in the form of the 

— * 

direction of the p integral and the azimuthal angle of A; if p is taken as the z axis. 

— * 

The dot product p-k—pk cos 9pk and 

g2 ^ = p2 ^ /j2 ^ 2pk COS 9pk COS Opk = ^"T ■ (B2) 

ApK 

If we change variables from p. A;, cos^p^ to p, k, q we should differentiate the above holding 
p, k fixed, giving 

q 

dcos 9pk — —dq. (B3) 

'pKi 
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Therefore we can rewrite the integration as 



d^pd?k 1 
(27r)6 " 2%4 



/•oo /-oo rp+k POD POO Mp+k)"^ 

/ pdp I kdk / gc^? = — — / dp^ / dk"^ I dq^ (B4) 

./o Jo J\p-k\ 2 TT Jo io J{p-k)^ 



which are a convenient set of integration variables. In particular, all dot products we will 
encounter can be written directly in terms of the integration variables using 

^ ^ k'^-p^-q^ - ^ p^-k'^-q^ ^ - q^ - p^ - k"^ 
-p-q= ^ , -k-q^ , p-k^ . (B5) 

The remaining two-loop diagram, the figure-8, can be performed using the same integra- 
tion variables; the two lines have momentum p and k, so the q integral can be done directly, 
J qdq = 2pk. This sort of reduction always works, because we can always consider a 4-point 
vertex to be two three-point vertices connected by a propagator, with the propagator col- 
lapsed to a point. So diagrams containing 4-point vertices can be written with the same 
variables as the diagram containing this "collapsed" propagator. 



2. Three loops: Mercedes 

Now we seek a similar set of integration variables for the Mercedes diagram. 




P 

Note that I = k + q and similarly k' — —p — k and — P — q- The phase space is 
determined by the triple integral 

fd'pd^kd^q r 2> ,2,, 2^ IStt' f . a f\ a f\j. ma^ 
J — (27r)9 — Pdpkdkqdqj^^ J dcosOpk J d cos 9pq J d(ppk-pq (B6) 

where is the azimuthal angle between the pk plane and the pq plane, we have used the 
symmetry of the (f) integration to reduce it from [0, 27r] to [0, tt], and Eulerian angles have 
again been performed. 

Using the same trick as before, we can rewrite this integral as 

-j^ POO POO POO rp+k rp+q 

/ dp kdk / qdq / k'dk' / qdq / d(t)pq-pk ■ (B7) 

^ Jo Jo Jo J\p~k\ J\p-q\ Jo 

and we would like to rewrite the integral as an I integral. To do so, write out an expression 
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for 



= (k + q)' = k' + q' + 2k- g , 

k ■ q = kq (cos 9pq cos 9pk + sin 9pq sin 9pk cos 0) = 



/CQ' cos 9pq COS 9pk 



{p'' + q^ - q'^){k'^ - - e) 
4p2 



/eg sin ^pq sin ^pfe = ^ {k'^ — k"^ cos^ 6pk) {q^ — q^ cos^ 6pq) 



and hence 



COS0 



{2p^^+2p^q'^+2q'^q''^-p^-q^-q''^){2p^k^+2p'^k'^+2k^k'^-p'^-k'^-k'^) 



(B8) 
(B9) 

(BIO) 
(Bll) 

(B12) 



Since the range of cos0 is from —1 to +1, we find that the range of at fixed p,k,q,k',q' 
is between 



1 

V 



{p\e+q^+k'^+q'') + q^k'^ + A;^^'' - p^ - k^ - k'\'') (B13) 

±"\/(V?H2p2g'2+2g2g'2_p4_^4_^/4^)^2p2A;2+2p2A;'2+2A;2A;'2-p4_fc4_^/4 



where the +(— ) sign represents the maximum (minimum) allowed value of P. 
Differentiating the expression for cos0 holding p, k, q, k', q' fixed, we find 



sin (pd4> — 



ApHdl 



where •>/ (...)(...) is the same long square root in the above expressions. Therefore 



(B14) 



(B15) 



sin(/>v/(::oR 

Writing sin = a/1 — cos2 and after significant algebra we find 

(B16) 

X = p2;2(^,.2^^,y2^^2^^/2_^2_^2)^^2^/2^^2^^/2^^2^^2_^2_^/2^ 

+eq''{q^+k''+p^+f-k-'-q'') - k^k'^-q^q'^-k^qH^-k'^'^P . (B17) 

Note that the expression for X has a symmetry, if hard to see. The momenta are in three 
pairs; {p,l), {q,k'), and {q',k) which are "opposite" momenta in the diagram (momenta 
which do not touch at a vertex). The first terms involve pairs of "opposite" momenta, the 
last terms involve triples of momenta meeting at a vertex. 
The total integration becomes 



pdp kdk qdq k'dk' q'dq' Idl 



(B18) 
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with integration limits listed previously. We have not written the integration limits in a 
symmetric way, but they are symmetric. 

The dot product of a pair of momenta which share a vertex are of form 

p-k= , p-q = (B19) 

where the sign difference is because in the first case the momenta are both directed out of 
the vertex while in the latter case one momentum enters and one exits the common vertex. 
For momenta which do not share a vertex, the dot product is 

p.l = p.{k + q} = ^—^ ^— (B20) 

and similarly for k-(f and q-k'. (For a mnemonic, note that q, k' are going from the beginning 
of one line to the end of the other; k, q' connect the beginnings of each line or the ends of 
each line.) We see that all dot products, including those for momenta on lines which do not 
meet at a vertex, have simple expressions in terms of momenta on lines. 

As mentioned before, we can use the same integration variables for 3-loop diagrams with 
one or more 4-point vertices. For instance, when the / propagator is collapsed into a 4-point 
vertex, one can immediately do the I integral; 

" (B21) 



X 2p 

However, if the integrand involves dot products which depend on / then we cannot do the 
/ integral immediately; we should instead interpret it as an angular integration which does 
not change the magnitudes of any momenta on the remaining lines, but which does affect 
some of their dot products. 
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TABLE II: ^ = 0.0 variational coefficients 
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